library(LinQseq)

# load UMI
umi.table <- read.table("/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895282_Hashtag-RNA.umi.txt.gz", 
    header = TRUE, row.names = 1)
umi.table <- as.sparse(umi.table)
umi.table <- CollapseSpeciesExpressionMatrix(object = umi.table)
# load HTO
hto.table <- read.csv(file = "/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895283_Hashtag-HTO-count.csv", 
    header = TRUE, row.names = 1)
remove.index <- which(rownames(hto.table) %in% c("bad_struct", "no_match", "total_reads"))
hto.table <- hto.table[-remove.index, ]

# load ADT
adt.table <- read.csv(file = "/Users/leil/Documents/Projects/scRNAseq/Multi-Modal/ECCITEseq/GSE108313_RAW/GSM2895284_Hashtag-ADT2-count.csv", 
    header = TRUE, row.names = 1)
remove.index <- which(rownames(adt.table) %in% c("bad_struct", "no_match", "total_reads"))
adt.table <- adt.table[-remove.index, ]

# only keep cells with RNA, ADT and HTO
joint.bcs <- intersect(colnames(umi.table), colnames(hto.table))
joint.bcs <- intersect(colnames(adt.table), joint.bcs)

umi.table <- umi.table[, joint.bcs]
hto.table <- as.matrix(hto.table[, joint.bcs])
adt.table <- as.matrix(adt.table[, joint.bcs])

rownames(x = hto.table) <- gsub(pattern = "-.+", replacement = "", x = rownames(x = hto.table))
rownames(x = adt.table) <- gsub(pattern = "-.+", replacement = "", x = rownames(x = adt.table))

pbmc <- CreateSeuratObject(counts = umi.table)

# add HTO
pbmc[["HTO"]] <- CreateAssayObject(counts = hto.table)
# Normalize HTO data, here we use centered log-ratio (CLR) transformation

# add ADT
pbmc[["ADT"]] <- CreateAssayObject(counts = adt.table)
# Normalize ADT data, here we use centered log-ratio (CLR) transformation

pbmc <- dataNormalization(object = pbmc)
pbmc <- HTODemux(pbmc, assay = "HTO", positive.quantile = 0.99)

table(pbmc$HTO_classification.global)
Idents(pbmc) <- "HTO_maxID"
VlnPlot(pbmc, assay = "HTO", features = rownames(pbmc[["HTO"]]), ncol = 3, pt.size = 0.01)

Idents(pbmc) <- "HTO_classification.global"

# t-SNE visualization
pbmc <- subset(pbmc, idents = "Negative", invert = TRUE)
# Calculate a distance matrix using HTO
hto.dist.mtx <- as.matrix(dist(t(GetAssayData(object = pbmc, assay = "HTO"))))

# Calculate tSNE embeddings with a distance matrix
pbmc <- RunTSNE(pbmc, distance.matrix = hto.dist.mtx, perplexity = 100)
DimPlot(pbmc)

DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne")

HTOHeatmap(pbmc, assay = "HTO", ncells = 5000)

pbmc <- FindVariableFeatures(object = pbmc, selection.method = "mean.var.plot")  # directly use Seurat Function
pbmc <- dataScaling(object = pbmc)

# PCA
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), verbose = FALSE)  # directly use Seurat Function

pbmc <- jointDistance(object = pbmc)
pbmc <- tsneFromDistane(object = pbmc, assay = "All")
pbmc <- clusteringFromDistance(object = pbmc, assay = "All")

# save(pbmc,file = 'HTO.rData')

plots <- generateGridDimPlot(pbmc, legend = FALSE, reduction.prefix = "tsne_", 
    darkTheme = FALSE)
listPlot(object = plots)

# Heatmap for joint clusters
heatMapPlot(object = pbmc, group.by = "jointClusterID")

# Heatmap for RNA clusters
heatMapPlot(object = pbmc, group.by = "rnaClusterID")

# Heatmap for ADT clusters
heatMapPlot(object = pbmc, group.by = "adtClusterID")

phto1 <- DimPlot(pbmc, group.by = "HTO_classification.global", reduction = "tsne_joint") + 
    DarkTheme()
phto2 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_joint") + DarkTheme()
phto3 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_rna") + DarkTheme()
phto4 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_adt") + DarkTheme()
plot_grid(phto1, phto2, phto3, phto4, ncol = 2)

phto1 <- DimPlot(pbmc, group.by = "HTO_classification.global", reduction = "tsne_rna") + 
    xlab("t-SNE1") + ylab("t-SNE2")
phto2 <- DimPlot(pbmc, group.by = "HTO_maxID", reduction = "tsne_rna") + xlab("t-SNE1") + 
    ylab("t-SNE2")
plot_grid(phto1, phto2, ncol = 2)

features.adt <- rownames(pbmc@assays[["ADT"]])
features.hto <- rownames(pbmc@assays[["HTO"]])

FeaturePlot(pbmc, features = features.adt, min.cutoff = "q05", max.cutoff = "q95", 
    ncol = 3, reduction = "tsne_adt")

FeaturePlot(pbmc, features = features.hto, min.cutoff = "q05", max.cutoff = "q95", 
    ncol = 3, reduction = "tsne_adt")

FeaturePlot(pbmc, features = features.adt, min.cutoff = "q05", max.cutoff = "q95", 
    ncol = 3, reduction = "tsne_rna")

FeaturePlot(pbmc, features = features.hto, min.cutoff = "q05", max.cutoff = "q95", 
    ncol = 3, reduction = "tsne_rna")